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Abstract. The electronic ground state of a periodic crystalline solid is usually de- 
scribed in terms of extended Bloch orbitals; localized Wannier functions can alterna- 
tively be used. These two representations are connected by families of unitary trans- 
formations, carrying a large degree of arbitrariness. We have developed a localization 
algorithm that allows one to iteratively transform the extended Bloch orbitals of a 
first-principles calculation into a unique set of maximally localized Wannier functions. 
We apply this formalism here to the case of cubic BaTiC>3. The purpose is twofold. 
First, a localized-orbital picture allows a meaningful band- by-band decomposition of 
the whole Bloch band complex. In perovskites, these Wannier functions are centered 
on the atomic sites and display clearly a s, p, d, or hybrid character. Second, since 
the centers of the Wannier functions map the polarization field onto localized point 
charges, the ground state dielectric properties become readily available. We study the 
Born effective charges of the paraelectric phase of BaTiC>3. We are able to identify not 
only the contributions that come from a given group of bands, but also the individual 
contributions from the "atomic" Wannier functions that comprise each of these groups. 



INTRODUCTION 

The electronic ground state of a periodic solid, in the independent particle ap- 
proximation, is naturally labelled according to the prescriptions of Bloch's theorem: 
single-particle orbitals are assigned a quantum number k for the crystal momentum, 
together with a band index n. Although this choice is widely used in electronic 
structure calculations, alternative representations are available. The Wannier rep- 
resentation [1] , essentially a real-space picture of localized orbitals, assigns as quan- 
tum numbers the lattice vector R of the cell where the orbital is localized, together 
with the band index n. Wannier functions can be a powerful tool in the study of the 
electronic and dielectric properties of materials: they are the solid-state equivalent 
of "localized molecular orbitals" [2], and thus provide an insightful picture of the 
nature of chemical bonding, otherwise missing from the Bloch picture of extended 
orbitals. In addition, the modern theory of polarization [3] directly relates the 
centers of the Wannier functions to the macroscopic polarization of a crystalline 
insulator. 



Wannier functions are strongly non-unique. This is a consequence of the phase 
indeterminacy e 1 ^ 71 ^ that Bloch orbitals t/j n k have at every wavevector k. This 
indeterminacy is actually more general than just the phase factors: Bloch orbitals 
belonging to a composite group of bands (i.e. bands that are connected between 
themselves by degeneracies, but separated from others by energy gaps) can un- 
dergo arbitrary unitary transformations between themselves at every k. We 
have recently developed a procedure [4] that can iteratively refine these otherwise 
arbitrary degrees of freedom, so that they lead to Wannier functions that are well 
defined and that are localized around their centers (in particular, they minimize 
the second moment around the centers). Such a procedure can be applied either 
to a whole band complex of Bloch orbitals, or just to some isolated subgroups. 

As a natural first application of this technique, we present here results for the 
case of BaTiC>3 in the cubic phase. Perovskite ferroelectrics, of which BaTiC>3 is a 
paradigmatic example, owe their very rich phenomenology to the subtle competi- 
tion of several degrees of freedom, balancing the long-range dipole-dipole interac- 
tion with short-range Pauli repulsion. One of the striking features is the display 
of anomalously large Born effective charges [5]. Their origin is understood in a 
simple tight-binding picture [6]: the change in the bond length (Ti-0 in this case) 
corresponds to a dynamic charge transfer that is stronger when the bonding is bor- 
derline between ionic and covalent. Localized Wannier functions can thus be used 
to investigate the nature of this bonding, and to monitor the changes that follow 
a ferroelectric distortion. Additionally, the displacement of each Wannier center 
relates directly to the effective charge contribution of its orbital, and can be used to 
identify the nominal and the anomalous contributions to the polarization induced 
by an atomic displacement. 



Electronic structure calculations are carried out using periodic boundary condi- 
tions. This is the most natural choice to study perfect crystals and to minimize 
finite size-effects in the study of several non-periodic systems (e.g. surfaces, or 
impurities). The one-particle effective Hamiltonian H then commutes with the 
lattice-translation operator Tr, allowing one to choose as common eigenstates the 
Bloch orbitals | i/j n k ), 



where « n k(r) has the periodicity of the Hamiltonian. There is an arbitrary phase 
n (k), periodic in reciprocal space, that is not assigned by the Schrodinger equation 
and that we have written out explicitly. We obtain a (non-unique) Wannier repre- 
sentation using any unitary transformation of the form ( nk | Rn ) = e lipn ^ _dc R : 



METHOD 




(1) 





Here V is the real-space primitive cell volume. It is easily shown that the | Rn ) form 
an ort honor mal set, and that two Wannier functions | Rn) and | R'n) transform 
into each other with a translation of a lattice vector R — R' [7] . The arbitrariness 
that is present in ip n (k) [or n (k)] propagates to the resulting Wannier functions, 
making the Wannier representation non-unique. Since the electronic energy func- 
tional in an insulator is also invariant with respect to a unitary transformation of its 
n occupied Bloch orbitals, there is additional freedom associated with the choice 
of a full unitary matrix (and not just a diagonal one) transforming the orbitals 
between themselves at every wavevector k. Thus, the most general operation that 
transforms the Bloch orbitals into Wannier functions is given by 

\ Rn ) = j BZ E UgH I ) e~ lk - R dk . (3) 

The Wannier functions w n (r — R) = | Rn), for non-pathological choices of phases, 
are "localized": for a R« far away from R, w n (TLi — R) is a combination of terms 
like J BZ -u m k(0)e lk '( Ri ~ R ) dk, which are small due to the rapidly varying character 
of the exponential factor [7] . 



Maximally-localized Wannier functions 

Several heuristic approaches have been developed that construct reasonable sets 
of Wannier functions, reducing the arbitrariness in the with symmetry consid- 
erations and analiticity requirements [8], or explicitly employing projection tech- 
niques on the occupied subspace spanned by the Bloch orbitals [9]. At variance 
with those approaches, we introduce a well-defined localization criterion, choosing 
the functional 

(4) 

n 

as the measure of the spread of the Wannier functions. The sum runs over the n 
functions | On); (r 2 ) n and f n = (r) n are the expectation values (On | r 2 | On) and 
( On | r | On ). Given a set of Bloch orbitals | ip m k ), the goal is to find the choice of 
U^l in (3) that minimizes the values of the localization functional (4). We are able 
to express the gradient G = of the localization functional with respect to an 
infinitesimal unitary rotation of our set of Bloch orbitals 

Kk) -> kO + E^iKk) , (5) 

m 

where dW an infinitesimal antiunitary matrix dW' = —dW such that 

77( k ) = 8 + dW^ (6) 

This provides an equation of motion for the evolution of the UjV, and of the | Rn ) 
derived in (3), towards the minimum of Q; small finite steps in the direction opposite 
to the gradient decrease the value of Q, until a minimum is reached. 



1 Real-space representation 



There are several interesting consequences stemming from the choice of (4) as the 
localization functional, that we briefly summarize here. Adding and subtracting 

~ 2 

the off-diagonal components Q = J2 n ^2nm^on (R?7i|r|0n) , we obtain the decom- 
position Q = Qi + fin + floB, where Qi, Qd and Qqb are respectively 



(r 2 ) n -£|(Rm|r|On> 

Rm 



fi D = EE \(Rn\r\On) 

n R^O 

^od = E E|( Rm l r l° n ) 

m^n R 

It can be shown that all terms are positive- definite (in particular Qi, see Ref. [4]); 
more importantly, fl \ is also gauge-invariant, i.e., it is invariant under any arbitrary 
unitary transformation (3) of the Bloch orbitals. The minimization procedure thus 
corresponds to the minimization of Q = Qd + Ood- At the minimum, the elements 

2 

(Rm|r|On) are as small as possible, realizing the best compromise in the simulta- 
neous diagonalization, within the space of the Bloch bands considered, of the three 
position operators x, y and z (which do not in general commute when projected 
within this space). 



2 Reciprocal- space representation 

As shown by Blount [7], matrix elements of the position operator between Wan- 
nier functions take the form 

(Rra|r|Om) = i j dke lkR (« nk | V k \u mk ) (7) 

and 

(Rn\r 2 \0m) = J dk e^ (u nk \V 2 k \u mk ) . (8) 

These expressions provide the needed connection with our underlying Bloch for- 
malism, since they allow us to express the localization functional Q in terms of the 
matrix elements of Vk and V k . We thus determine the Bloch orbitals \u mk ) on a 
regular mesh of k-points, and use finite differences to evaluate the above deriva- 
tives. For any given k-point in a regular cubic mesh (sc, fee, bec), we have a star 
b of Z k-points that are first-neighbors; their weights in the evaluation of deriva- 
tives are Wb = 3/Zb 2 . We define M^ b ) = (u mk \u nik+ b) as the matrix elements 



between Bloch orbitals at neighboring k-points. The M^ b ) are a central quantity 
in our formalism, since we can then express all the contributions to the localization 
functional using the connection made by Blount, together with our finite-difference 
evaluations of the gradients. After some algebra we obtain [4] 

fi I=^E^ f^bands " E l^i k n b) | 2 ) , (9) 
JV k,b \ run ) 

n D = ^E^E l M ^ b) l 2 , (io) 

k,b rriy^n 

and 

fi D = ^E^E (" Im ^ M in h) ~ Im \nMt b) ) 2 . (11) 

^ V k,b n ^ ' 

From these, we can calculate the change in the localization functional in response 
to an infinitesimal unitary transformation of the Bloch orbitals, as a function of the 
M^'; once these steepest-descents are available, it is straightforward to construct 
a procedure that updates the towards the minimum of the functional. 

RESULTS: THE CASE OF BATIOg 

We study here the cubic phase of BaTi0 3 , using a plane- wave total-energy 
pseudopotential approach with the local-density approximation to the exchange- 
correlation functional. We use norm- conserving pseudopotentials in the Kleinman- 
Bylander representation, with q c kinetic-energy tuning [10] for the oxygen atom and 
a Troullier-Martins procedure [11] for the titanium, to bring the cutoff convergence 
down to 900 eV. The 3s and 3p levels of titanium have been included in the valence. 
The Brillouin zone is sampled with a 4 x 4 x 4 Monkhorst-Pack mesh; the lattice 
parameter used is 3.98A. 

Wannier functions of cubic BaTiC>3 

The minimization of the total energy provides the Kohn-Sham Bloch orbitals on 
a regular mesh of k-points, that are then used as a starting point for the construc- 
tion of the Wannier functions. The subsequent minimization of the localization 
functional determines the that correspond to the maximally-localized Wan- 
nier functions. In BaTi0 3 there are several groups of bands that are separated by 
gaps. In order of increasing energy, we have the band groups corresponding to the 
Ti 3s (1), Ti 3p (3), Ba 5s (1), O 2s (3), Ba hp (3) and O 2p (9) levels (in paren- 
thesis are the number of bands in each group). We initially consider each group of 



FIGURE 1. Left panel: oxygen-centered Wannier function from the O 2s 3-band group (the 

atom is surrounded by four Ba atoms on the sides, and two Ti atoms on top and bottom). 
Center and right panels: barium-centered Wannier function from the Ba 5p 3-band group, and 
barium-centered Wannier function from the Ba 5p and the O 2s 6-band group (the Ba atom is 
surrounded by 12 oxygens). 

bands separately, and perform the minimization on the 6 subspaces of dimensions 

1 x 1, 3 x 3, 1 x 1, 3 x 3, 3 x 3, and 9 x 9. 

The Wannier functions determined from the the Ti 3s or the Ti 3p groups strongly 
resemble atomic orbitals, slightly deformed by the underlying crystal potential, and 
are not shown here. We show instead in the left panel of Fig. 1 one of the three 
oxygen-centered Wannier functions that are derived from the oxygen 2s bands. In 
each unit cell there are three such functions, sitting on each of the oxygens. The 
orbital shows its atomic s character; there are some contributions on the Ti, with p 
and/or d character (the titanium is slightly embedded in a d z i orbital). Interesting 
results emerge from the localization of the Ba 5p bands (Fig. 1, center and right 
panels). These bands correspond to three orbitals in each unit cell, all centered 
on the barium and oriented along the three cristallographic directions (p x , p y and 
p z ). It can easily be seen (center panel) that, in addition to the distinctive atomic 
p orbital on the barium, there are significant sp-like contributions sitting on 8 
of the 12 neighboring oxygens. This supports the suggestion that barium in this 
compound has some covalent character [12], and it is consistent with the anomalous 
effective-charge contributions that come from this group of Wannier functions (see 
next subsection). It is interesting to note that if we decide to treat the Ba hp bands 
together with the O 2s bands, we can (obviously) increase the degree of localization 
of each orbital. In this latter case (right panel) the sp contributions on the oxygens 
decrease, being transferred to the 2s orbitals localized on the oxygens themselves. 

Finally, we examine the 9 oxygen 2p bands that result from the hybridization of 
the O 2p electrons with the Ba 6s and the Ti 4s and 3d. There are three localized 
orbitals on each oxygen, oriented along the Ti-O-Ti bonds. We label two of these 
orbitals as 7r and one as a, according to their symmetry along the Ti-O-Ti axis. One 
of the a orbitals is shown in Fig. 2, first with the atoms in their ideal positions (left 




FIGURE 2. Oxygen-centered a Wannier function from the localization of the O 2p 9-band 
group. The orbital is oriented along the Ti-O-Ti bond; the four Ba atoms neighboring the central 
oxygen and the two Ti atoms on top and bottom are also shown. Left panel: ideal atomic 
positions. Right panel: same, but with the titanium atoms displaced downwards. 



panel) and then with the Ti atoms displaced along the Ti-0 bond (right panel). 
The a orbital shows clearly the hybridization between the oxygen p orbital oriented 
along the Ti-O-Ti direction and the d z i orbital of the titanium. The a and even 
more the 7r orbitals have strong anomalous contributions to the effective charges, 
that can be visualized with the large charge transfer from the oxygen atoms in 
response to the titanium displacements. 

Band-by-band decomposition of the Born effective charges 

The Born dynamical effective charges describe the change in macroscopic polar- 
ization that is induced by the displacement of a given ion. As such, they play a 
fundamental role in determining the dynamical properties of insulating crystals, 
and are a powerful tool to investigate the dielectric and ferroelectric properties of 
materials. They also determine the splitting of the infrared-active optical modes; 
in simpler compounds (e.g. GaAs) they can be unambiguously determined from 
the experimental phonon dispersions. Perovskite ferroelectrics display anomalous, 
large effective charges, that can be almost double their nominal ionic value. The 
origin of this effect lies in the large dynamical charge transfer that takes place when 
moving away from the high-symmetry cubic phase (i.e., going from more ionic to 
more covalent bonding). Orbital hybridization is necessary for this transfer to take 
place; for this reason, our localized-orbitals approach provides an insightful tool in 
examining these effects. In the language of the modern thory of polarization [3] , the 
anomalous contribution is determined by the relative displacement of the Wannier 
centers with respect to the ion that is being moved. If the bonding were purely 
ionic, electrons (and thus Wannier centers) would be firmly localized on each anion, 
and move rigidly with it. This is not the case in perovskite oxides. The anomalous 
contribution is often traced [13] to substantial hybridization between the oxygen p 
orbitals and the d orbitals of the atom in the B site (Ti, in this case). The picture 
can be somewhat more complex, with other group of bands playing a role in the 
anomalous dielectric behavior. 

We present in Table 1 a full decomposition of the effective charges in BaTi0 3 
coming from the different groups of bands; and, inside each group, coming from the 
individual Wannier functions identified by the localization procedure. We compare 
the results for the groups of bands with those obtained by Ghosez et al. (Ref. [12]). 
We find very good agreement, given the difference in the pseudopotentials used 
and our choice of lattice parameter. These results underline the conclusion that 
the decomposition into band groups is consistently defined in the linear-response 
formalism only when the calculations are performed in the so-called diagonal gauge 
[12]. The effective-charge tensor reduces to a scalar for barium and titanium, while 
it is diagonal with two inequivalent components (those parallel and perpendicular 
to the Ti-0 bonds) for the oxygens. Several facts stand out from an inspection 
of the table. The anomalous effective charges originate not only from the oxygen 
2p bands; there are also sizeable contributions originating from the barium hp 



TABLE 1. Born effective charges decomposed into atomic contributions, for BaTiC>3. The Z* ; are calculated 
by displacing a Ba, Ti or a O sublattice (Oi is parallel to the Ti-0 direction, O2 perpendicular) by Sa , where 
S = 0.002. The numbers in parenthesis are those of Ghosez et al. (ao=3.94A). Horizontal lines group band 
complexes that have been treated together in the Wannier minimization. Top line in each group is the total 
for that group. We have added at the bottom the O 2p decomposition obtained in a calculation performed 
using a Ti atom with the 3s and 3p levels frozen in the core (Ti 4+ ). 
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orbitals and even from the oxygen 2s. More notably, there is a wide range of 
compensating effects between groups of bands and between different orbitals inside 
each group. The partial cancellation of such large orbital polarizabilities hints again 
at the complexity of these materials, which exhibit such a wide range of equilibrium 
properties as a consequence of the existence of many of these competing effects. 
Titanium shows the strongest deviations from a naive ionic picture. The a and 
7i oxygen orbitals carry a positive electronic Z* contribution, equal respectively to 
0.81 and 1.78. It should be noted that it is always the O 2p ix orbitals that carry 
the largest anomalous charge. The O 2p contributions to the Oi effective charges 
(i.e. in the direction of the Ti-0 bond) are also anomalous, up to —1.10 for each n 
orbital (in addition to the nominal —2.00 for each orbital). 

Finally, we present at the bottom of Table 1 the oxygen 2p decomposition per- 
formed in a calculation where the 3s and 3p orbitals of the titanium have been 
removed from the pseudopotential (i.e., removed from the valence and frozen in the 
core). It is interesting to note that most contributions are completely unchanged; 
differences arise only for the anomalous contributions for the Ti and the Oi dis- 
placements (whose anomality, incidentally, is reduced if this more approximate 
formalism is employed). 

CONCLUSIONS 

We have summarized here our formalism for obtaining maximally- localized Wan- 
nier functions from the Bloch orbitals of an ab-initio electronic structure calcula- 
tion. This formalism can be very helpful in understanding the chemical and dielec- 
tric properties of materials. Perovskite ferroelectrics are a particularly promising 
class of systems to be studied, since the nature of the bonding and hybridization 
can have a striking influence on the dielectric properties and on the development 
of ferroelectricity. At variance with other approaches, our method allows for a 
decomposition of electronic properties (e.g., the effective charges) into meaningful 
atomic contributions. In the case of BaTiC>3, it elucidates in particular the origins 
of the large anomalous contributions to the effective charges. 
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